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Conventional methods to calculate the thermodynamics of crystals evaluate the harmonic phonon 
spectra and therefore do not work in frequent and important situations where the crystal structure 
is unstable in the harmonic approximation, such as the body-centered cubic (bcc) crystal structure 
when it appears as a high-temperature phase of many metals. A method for calculating temperature 
dependent phonon spectra self consistently from first principles has been developed to address this 
issue. The method combines concepts from Born's inter-atomic self-consistent phonon approach 
with first principles calculations of accurate inter-atomic forces in a super-cell. The method has 
been tested on the high temperature bcc phase of Ti, Zr and Hf, as representative examples, and is 
found to reproduce the observed high temperature phonon frequencies with good accuracy. 



Many elements, alloys, and compounds appear in crys- 
tal structures which should not be energetically stable. 
The inter-atomic interactions places these systems at en- 
ergy saddle points on the potential surface for atomic 
positions corresponding to the lattice sites of these struc- 
tures rather than minima for statically stable structures. 
The body centered cubic (bcc) structure prevails as the 
simplest and best known example. Although a stable 
structure at low temperatures for several elements in the 
Periodic Table, bcc becomes unstable in the harmonic 
approximation [H, [3, Q for the group IVB elements, for 
the rare-earth elements, for the actinides, and for several 
alkaline-earth elements. Nevertheless, at elevated tem- 
peratures the bcc structure emerges as the stable crystal 
structure for all these elements. Zener considered this 
enigma long ago and proposed a possible explanation: 
the large vibrational entropy of the bcc crystal structure 
makes it thermodynamically favourable at finite temper- 
atures [H. Also, Grimvall et al 0] pointed out the impor- 
tance of electronic entropy in the stabilization of the bcc 
crystal structure of the group IVB elements Ti and Zr. 

So far no satisfactory, quantitative explanation has 
been presented for this situation. Density functional the- 
ory (DFT) [|| forms the basis of contemporary micro- 
scopic solid state theory and allows, in principle, to calcu- 
late different properties of crystals completely ab initio, 
without any fitting parameters. In particular, phonon 
spectra in the harmonic approximation can be efficiently 
evaluated in this way However, for the bcc phases 
mentioned above the phonon spectra in the harmonic ap- 
proximation reveal imaginary phonon frequencies of e.g. 
Zr[|[i for some wave- vectors, which shows that the bcc 
phase is from a lattice dynamics point of view unstable 
(hence these elements are energetically unstable, and are 
referred to as dynamically unstable in the bcc phase). 
A straightforward calculation using DFT molecular dy- 
namics (MD) should in principle be able to reproduce 
the stability of the bcc phase for the above discussed el- 



ements, since MD implicitly include temperature effects. 
However, MD suffers from the fact that obtaining reliable 
free energies implies a computationally very demanding 
task, which in many cases make these types of calcula- 
tions intractable. 

We propose here a solution to this problem, which 
builds on a self-consistent ab initio lattice dynamics 
(SCAILD) approach. In this paper we describe the es- 
sential aspects of our method and apply it to the problem 
of stability of the bcc phase for the group IVB elements. 
Although several aspects of our proposed theory have not 
been considered before we note that it conceptually has 
similarities with the self-consistent phonon approach by 
Bor n II ill , and that several other self-consistent methods 
lil ljjhave been developed in the past. We will show 
that the SCAILD theory gives phonon spectra of the bcc 
phase of Ti, Zr and Hf which are in agreement with ob- 
servations 15, 0, 17 1 . Although we will in the rest of 
this manuscript focus on the group IV elements, we point 
out here that what we provide is a general scheme which 
can be used for any element and compound. 

Self consistent phonon calculations are a natural exten- 
sion of the theory of the harmonic lattice, and we initi- 
ate our methodological description by first presenting the 
most important features of this theory. The Hamiltonian 
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describes a harmonic lattice where R are the equilib- 
rium lattice positions of the atoms, Ur the displacements 
of the atoms, PR_the momentum of the atoms, M the 
atomic mass and <E> the inter-atomic force constant ma- 
trices (Here the vectors R refere to the positions of a 
Bravais lattice). Diagonalizing the dynamical matrix 
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for each wave vector k in the first Brillouin zone one 
finds the eigenvalues Wks and eigenvectors eks of different 
phonon modes (longitudinal or transverse) labeled by the 
symbol s, N being the number of atoms. Introducing the 
canonical phonon coordinates Ur and Pr 



ks^kse 



ikR 



ikR 



(3) 
(4) 



k,s 



allows a separation of the original Hamiltonian of the 
crystal into the Hamiltonians of 3N independent har- 
monic oscillators. 

The thermodynamic average of the operators Q^Qks 
determines the mean-square atomic displacements and is 
given by 



(QlQk S ) = -f- 
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where n(x) — l/(e x — 1) is the Planck function. In the 
classical limit, i.e for sufficiently high temperatures, the 
operators (l/\/A7) Qts are replaced by real numbers, 



Aws = ±1 
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Calculating the gradient of the potential energy in 
Eqn. [1] with respect to the atomic displacements gives 
the restoring force 



F R = -^i>(R-R')U R , 



(7) 



Fourier transforming Eqn. [7] and substituting Ur with 
the expression in Eqn. 3 gives 
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Finally, using the orthogonality of the eigenvectors 6k s 
the phonon frequencies can be expressed as 
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The equations discussed so far can be solved for dy- 
namically stable materials, where each atom is located 
in a minimum of the function C/r. It is important to 
note that this does not have to correspond to a global 
total energy minimum of the lattice, a local minimum 
suffices. For dynamically unstable materials [/r does 
not have a minimum at the lattice sites of the crystal 
structure. In this situation the equations discussed so 
far can not be used straight forwardly since they result 
in imaginary phonon frequencies. This represents a situa- 
tion where the lattice under consideration spontaneously 



shifts atomic planes and/or atomic positions so that a 
new crystal structure lowers the total energy. We demon- 
strate the problem at hand by comparing in FigJT] the 
calculated (zero temperature) phonon spectra of the bec 
phase of the group IVB elements with experimental data 
obtained at elevated temperatures. At these tempera- 
tures the group IVB elements are observed to be stable 
in the bec crystal structure, and the measured phonon 
frequencies are naturally positive for all lattice vectors. 
Fig. 1 shows that calculations using a static bec lattice re- 
sult in a dynamically unstable situation with imaginary 
phonon frequencies. It should be noted that the fail- 
ure describing the bec phase of the group IVB elements 
using harmonic lattice theory (Fig.l, right column) is 
not caused by any obvious error in the energy functional 
used, and are likely not to be improved even if an exact 
functional for a static lattice were found. 

In order to properly describe the high temperature 
phase of the group IVB elements on must include the 
interaction between phonons [HI]. As a result, phonon 
frequencies turn out to be temperature dependent which 
we explore numerically in this study. However, we neglect 
the phonon damping due to decay processes of phonons 
(see, e.g., Ref. [19J and Refs. therein), another anhar- 
monic effect. In the present calculations thermal expan- 
sion effects have not been taken into account, all calcula- 
tions have been performed at constant volume. Further- 
more, the thermal excitations of the electronic subsystem 
has not been considered in the present calculations of the 
phonon frequencies. 

The method used to calculate temperature dependent 
phonons presented in this paper considers a supercell con- 
taining a number of atoms which are allowed to deviate 
from the lattice positions stipulated by the crystal struc- 
ture. The deviations are calculated as a function of tem- 
perature, by solving equations 3 to 9 self consistently. 
The deviation of the atomic positions away from the ideal 
lattice points provides an extra entropy to the system and 
the stabilization of the bec structure for the group IVB 
elements as a function of increasing temperature may as 
we will see below be found. 

As regards the calculational details of the force calcu- 
lation we used the VASP package within the gen- 
eralized gradient approximation (GGA). The PAW po- 
tentials used required energy cutoffs of 197 eV for Ti, 
175 eV for Zr, and 243 eV for Hf. The k-point mesh 
was a 6x6x6 Monkhorst-Pack grid, and the supercell used 
was obtained by increasing the bec primitive cell 4 times 
along the 3 primitive lattice vectors. 

In practice our calculations are done by first calcu- 
lating a starting guess for the phonon dispersions by 
means of a standard supercell calculation, see e.g Ref. 
I20I . The phonon frequencies corresponding to k-vectors 
commensurate with the supercell are then used to calcu- 
late the atomic displacements through Eqns. 3-6. Here it 
should be noted that the signs of the amplitudes Au s (see 
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FIG. 1: The phonon dispersions of the group IVB metals. The solid lines are the first principles self consistent phonon 
calculations. In the left column the finite temperature calculations, and in the right column the T = K calculations. The 
filled circles are the experimental data of Ref. [15l. [lril. . 
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FIG. 2: The change in free energy between two consecutive 
iterations, here plotted as a function of the number of itera- 
tions. The inset in the figure shows the same plots but at a 
smaller energy scale. 



introduce any extra approximation, it merely guarantees 
that the longitudinal and transverse modes are fixed to 
the modes of the bee lattice. Starting from the equilib- 
rium geometry used in the initial supercell calculation, 
the atoms are displaced according to Eq. 3, and the 
forces on these displaced atoms are calculated. From the 
Fourier transform of the atomic forces a new set of fre- 
quencies are calculated through Eqn ©. To retain the 
correct symmetry of the calculated phonon dispersion the 
symmetries of the different k- vectors are restored by 
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where 5(k) is the symmetry group of the wave vector k, 
and mk the number of elements of the group. From the 
different iterations frequency distributions of the modes 
are obtained, and a new set of frequencies are supplied 
by the mean frequencies of these distributions, 
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Eqnj6]), should be chosen randomly, with equal probabil- 
ities for + and -. This is an approximation to the pro- 
cedure in which A^s is sampled continuously to obtain 
the correct mean square deviations of the modes. This 
approximation however conserves correctly the property 
that the displacements in Eqn. 3 are R-dependcnt. Fur- 
thermore, it should be noted that the eigenvectors ek s cal- 
culated in the initial calculation are not updated through- 
out the rest of the procedure. This however does not 



where f2ks(*)i i = 1, -/V are the symmetry restored fre- 
quencies from all iterations. The new set of frequencies 
calculated in (jlip determine a new set of displacements 
used to calculate a new set of forces. Philosophically our 
approach is similar to Born's self consistent phonon the- 
ory, with the main difference being that we consider a 
direct force calculation from a super cell with Hellman- 
Feynman forces calculated from density functional the- 
ory. 
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Figure Q] shows the calculated phonon dis pers ions to- 
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gether with the experimental data of Ref. 
for the bcc phase of the group IVB metals at tempera- 
tures 1293 K, 1188 K, and 2073 K for Ti, Zr, and Hf, 
respectively. The finite temperature calculations predict 
the stability of the bcc phase of all group IVB metals 
by promoting the frequencies of the phonons along the T 
to N symmetry line and around the P symmetry point 
from imaginary to real. The finite temperature calcula- 
tions of phonons result in an overall quantitative agree- 
ment with experimental values. Smaller deviations are 
observed around the P and H point of the Brillouin-zone, 
most likely due to finite size effects of the supercell used 
in the calculations. 

From the self consistent phonon spectrum the free en- 
ergy is approximated from the density of states of the 
phonons g(w) through the expression 



approach has advantages over traditional methods such 
as MD simulations in that complications associated with 
metallic materials are avoided and, most importantly, 
that a much smaller set of atoms are needed. 
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which has been shown by Cochran et al. 2l| to give an 
entropy correct to leading order in anharmonic pertur- 
bation theory. Figure shows the convergence of free 
energy for the three elements considered in this work. In 
all calculations presented here the self consistent cycle 
was terminated when the difference in the approximate 
free energy of the lattice between two consecutive itera- 
tions was less than 1 meV. Convergence in the free energy 
with such accuracy is very encouraging and opens up the 
possibility to investigate temperature induced phase sta- 
bility for a very large set of materials, since the accuracy 
needed to e.g. resolve crystallographic energy differences 
is of the order of a few meV or more. This prediction has 
also been tested by using inter-atomic forces calculated 
with the embedded atom potentials of Ref. 2^, H- The 



free energy difference between the hep and bcc structures 
where calculated as functions of temperature for Ti and 
Zr. Here the theoretically predicted hep to bcc transition 
temperatures where within ~400 K of the corresponding 
experimental temperatures. 

In summary, a quantitative theory successfully ex- 
plains the long lasting question concerning thermal, en- 
tropy driven stabilization of dynamically unstable mate- 
rials. Application to the group IVB elements reproduces 
the measured phonon spectrum of these elements at el- 
evated temperatures with good accuracy. We note that 
the presented method reproduces observed high temper- 
ature phonon spectra with good accuracy and that the 
method when used at low, but non-zero temperatures, 
results in imaginary frequencies for e.g. bcc Ti. This 
shows that at low temperatures this element is unstable 
in the bcc phase, in agreement with observations. Other 
systems where one can expect success of this method are 
the bcc phase of f-electron materials as well as, the high 
pressure phase of Fe, and many of the ferroelectrics. The 
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